%matplotlib notebook
import numpy as np
import matplotlib as mpl
import mpl_toolkits.mplot3d
import matplotlib.pyplot as plt
import scipy.signal
import scipy.fftpack
mpl.rcParams["figure.dpi"] = 300
mpl.rcParams["figure.figsize"] = (10, 5)
import gc
gc.collect()
waveform0 = plt.imread("mouse1-processed.png")
waveform = np.dot(waveform0[..., :3], [0.299, 0.587, 0.114])[0: 644, 11: 1023]
waveform
np.max(waveform), np.min(waveform)
clear all
clc
tic
%图片读取
waveform0=rgb2gray(imread('Image1.bmp'));
waveform0=im2double(waveform0);
waveform0=waveform0(1:645,12:1024);
plt.imshow(waveform, cmap="gray")
plt.colorbar()
# plt.show()
waveform.shape
图片尺寸是644x1012
plt.imshow(waveform, cmap="jet")
plt.colorbar()
# plt.show()
fs = 300e+6
c0 = 1500
x0 = c0 / fs
t0 = 1 / fs
n = 30
nn = 2 * n + 1
fs=300e6; %采样频率200MHz 中心频率80MHz
c0=1500; %组织声速
x0=c0/fs; %每个采样点的距离
t0=1/fs; %采样间隔
n=30; %采样点数目2n+1
nn=2*n+1;
slope = np.zeros((waveform.shape[0] - 2 * n, waveform.shape[1]))
intercept = np.zeros((waveform.shape[0] - 2 * n, waveform.shape[1]))
f = list(range(0, n+1)) # f = [0, 1, 2, ..., n]
for i in range(0, waveform.shape[1] - 2):
print(i)
for j in range(0, waveform.shape[0] - nn):
temp = waveform[j: j + 2 * n + 1, i]
temp = temp * np.hamming(nn)
fft_temp = np.abs(np.fft.fft(temp))
FT = np.nan_to_num(20 * np.log10(fft_temp / np.max(fft_temp)))
# FT = fft_temp / np.max(fft_temp)
FT1 = FT[0: n + 1]
# print(FT1)
RE_a, RE_b = np.polyfit(f, FT1, 1)
slope[j, i] = RE_a
intercept[j, i] = RE_b
slope, intercept
slope=zeros(RF_Row-2*n,RF_Column);
intercept=zeros(RF_Row-2*n,RF_Column);
f=1:n+1;
for i=1:RF_Column-2
for j=1:RF_Row-nn
temp=waveform(j:j+2*n,i);
temp=temp.*hamming(nn);
fft_temp=abs(fft(temp));
FT=20*log10(fft_temp./max(max(fft_temp)));
FT1=FT(1:n+1);
[RE_a,RE_b]=line_reg1(FT1,f);
slope(j,i)=RE_a;
intercept(j,i)=RE_b;
end
end
toc
f = list(range(0, n+1))
def slopeFilter(X):
# temp = X * np.hamming(nn)
# fft_temp = np.abs(np.fft.fft(temp))
fft_temp = np.abs(np.fft.fft(X * np.hamming(nn)))
# fft_temp = np.abs(scipy.fftpack.fft(X * np.hamming(nn)))
FT = np.nan_to_num(20 * np.log10(fft_temp / np.max(fft_temp)))
# print(X)
return np.polyfit(f, FT[0: n+1], 1)[0]
# slope2 = scipy.ndimage.generic_filter(waveform, slopeFilter, size=(nn, 1))[nn-1:, 1-1:]
# slope2
%%timeit
slopeFilter(waveform[0: nn, 0])
slopeFilter需要约$400\,{\rm \mu s}$做一个窗口的计算。图片的大小是644x1012。所以总共有
$$
(644 - 61 + 1) \times 1012 = 591008
$$
个窗口。预计整张图片需要
$$
591008 \times 400 \cdot 10^{-6} = 236.4032 \,{\rm s}
$$
才能处理完。
希望能把单个窗口的运算时间压缩到$10 \,{\rm \mu s}$。这样每张图片需要的时间是 $$ 591008 \times 10 \cdot 10^{-6} = 5.91008 \,{\rm s} $$
%%timeit
for i in range(0, waveform.shape[0] - nn):
slopeFilter(waveform[i: i + nn, 0])
每列需要处理$644 - 61 + 1 = 584$个窗口。按照理论需要 $$ 584 \times 400 \cdot 10^{-6} = 0.2336 \,{\rm s} $$ 才能结束。然而实际测得需要$200 \,{\rm ms}$。
如果按照实际测得的时间算的话,整张图片需要 $$ 1012 \times 200 \cdot 10^{-3} = 202.4 \,{\rm s} $$ 来处理。
%%time
slope2 = np.zeros((waveform.shape[0] - nn + 1, waveform.shape[1]))
f = list(range(0, n+1))
for i in range(0, waveform.shape[1]): # column
# print(i)
for j in range(0, waveform.shape[0] - nn): # row
slope2[j, i] = slopeFilter(waveform[j: j + nn, i])
# slopeFilter(waveform[0: nn, 0])
slope2
# %time
用时$227\,{\rm s}$。
slopeFilter的延时¶%prun slopeFilter(waveform[0: nn, 0])
173 function calls in 0.001 seconds
Ordered by: internal time
ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 0.000 0.000 {built-in method numpy.fft.fftpack_lite.cfftf}
1 0.000 0.000 0.001 0.001 <ipython-input-11-84a0b36be0fd>:2(slopeFilter)
2 0.000 0.000 0.000 0.000 {built-in method numpy.linalg.lapack_lite.dgelsd}
1 0.000 0.000 0.000 0.000 linalg.py:1813(lstsq)
1 0.000 0.000 0.000 0.000 polynomial.py:398(polyfit)
1 0.000 0.000 0.000 0.000 function_base.py:3484(hamming)
13 0.000 0.000 0.000 0.000 {built-in method numpy.core.multiarray.array}
3 0.000 0.000 0.000 0.000 {method 'reduce' of 'numpy.ufunc' objects}
2 0.000 0.000 0.000 0.000 iostream.py:195(schedule)
发现np.fft比较慢,要$30 \,{\rm \mu s}$左右。如果用scipy.fftpack.fft代替的话,能快出$10 \,{\rm \mu s}$左右。
np.polyfit也是很慢的,要$180 \,{\rm \mu s}$左右。
numpy.fft改用scipy.fftpack.fft¶f = list(range(0, n+1))
def slopeFilter(X):
# temp = X * np.hamming(nn)
# fft_temp = np.abs(np.fft.fft(temp))
# fft_temp = np.abs(np.fft.fft(X * np.hamming(nn)))
fft_temp = np.abs(scipy.fftpack.fft(X * np.hamming(nn)))
FT = np.nan_to_num(20 * np.log10(fft_temp / np.max(fft_temp)))
# print(X)
return np.polyfit(f, FT[0: n+1], 1)[0]
# slope2 = scipy.ndimage.generic_filter(waveform, slopeFilter, size=(nn, 1))[nn-1:, 1-1:]
# slope2
所以这里把np.fft.fft替换成scipy.fftpack.fft。
%%timeit
slopeFilter(waveform[0: nn, 0])
%%timeit
for i in range(0, waveform.shape[0] - nn):
slopeFilter(waveform[i: i + nn, 0])
果然快了一点。
%%timeit
np.hamming(61)
每次生成Hamming窗都要花$16 \,{\rm \mu s}$。这部分时间可以省下来的。
f = list(range(0, n+1))
hammingWindow = np.hamming(nn)
def slopeFilter(X):
# temp = X * np.hamming(nn)
# fft_temp = np.abs(np.fft.fft(temp))
# fft_temp = np.abs(np.fft.fft(X * np.hamming(nn)))
fft_temp = np.abs(scipy.fftpack.fft(X * hammingWindow))
FT = np.nan_to_num(20 * np.log10(fft_temp / np.max(fft_temp)))
# print(X)
return np.polyfit(f, FT[0: n+1], 1)[0]
# slope2 = scipy.ndimage.generic_filter(waveform, slopeFilter, size=(nn, 1))[nn-1:, 1-1:]
# slope2
这里事先就生成一次Hamming窗,不再每次再在函数调用里生成一次。
%%timeit
slopeFilter(waveform[0: nn, 0])
%%timeit
for i in range(0, waveform.shape[0] - nn):
slopeFilter(waveform[i: i + nn, 0])
又省了$20 \,{\rm ms}$每列。
np.ALLOW_THREADS = 0
%%timeit
slopeFilter(waveform[0: nn, 0])
%%timeit
for i in range(0, waveform.shape[0] - nn):
slopeFilter(waveform[i: i + nn, 0])
没有什么作用。
import ipyparallel as ipp
rc = ipp.Client()
lview = rc.load_balanced_view()
# result.display_outputs()
lview
%%time
def task(i):
i, waveform = i
import numpy as np
import scipy.fftpack
fs = 300e+6
c0 = 1500
x0 = c0 / fs
t0 = 1 / fs
n = 30
nn = 2 * n + 1
# f = list(range(5, 29))
f = list(range(0, n+1))
hammingWindow = np.hamming(nn)
slope3 = np.zeros((waveform.shape[0] - nn + 1, 1))
intercept3 = np.zeros((waveform.shape[0] - nn + 1, 1))
def slopeFilter(X):
fft_temp = np.abs(scipy.fftpack.fft(X * hammingWindow))
FT = np.nan_to_num(20 * np.log10(fft_temp / np.max(fft_temp)))
# return np.polyfit(f, FT[5: 29], 1)
return np.polyfit(f, FT[0: n+1], 1)
for j in range(0, waveform.shape[0] - nn):
slope3[j], intercept3[j] = slopeFilter(waveform[j: j + nn])
return slope3, intercept3
# results = []
# f = list(range(0, n+1))
# for i in range(0, waveform.shape[1]):
# results.append(lview.apply(task, i, waveform[:, i]))
# print(i)
results = lview.map(task, [(i, waveform[:, i]) for i in range(0, waveform.shape[1])])
results
results.serial_time
results.wall_time
这是开4 workers的情况。用了$127 \,{\rm s}$。
results.serial_time
results.wall_time
这是开2 workers的情况。用了$166 \,{\rm s}$。
slope3 = np.hstack([i[0] for i in results.result()])
slope3
intercept3 = np.hstack([i[1] for i in results.result()])
intercept3
%time scipy.ndimage.generic_filter(waveform[0: 100, 0: 20], slopeFilter, size=(61, 1))
intercept
slope3.tofile("mouse1-slope.numpy")
intercept3.tofile("mouse1-intercept.numpy")
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
rows = np.array(range(slope2.shape[0]))
cols = np.array(range(slope2.shape[1]))
rows, cols = np.meshgrid(cols, rows)
ax.plot_surface(rows, cols, slope2, cmap="jet")
ax.set_aspect("equal")
# fig.colorbar(slope2)
slopeImage = mpl.imsh
plt.imshow(slope3, cmap="jet")
plt.colorbar()
# plt.show()
plt.imshow(20 * np.abs(slope3), cmap="jet")
plt.colorbar()
# plt.show()
plt.imshow(intercept3, cmap="jet")
plt.colorbar()
# plt.show()
plt.imshow(np.abs(intercept3), cmap="jet")
plt.colorbar()
# plt.show()
plt.imshow(slope3 * 16 + intercept3, cmap="jet")
plt.colorbar()
# plt.show()
colNormal = 600
rowNormal = 350
colTumor = 380
rowTumor = 320
rows = 80
cols = 100
rowOffset = 30
plt.imshow(waveform, cmap="jet")
plt.colorbar()
# plt.show()
plt.gca().add_patch(
plt.Rectangle([colNormal, rowNormal], cols, rows, fill=None, edgecolor="black") # normal
)
plt.gca().add_patch(
plt.Rectangle([colTumor, rowTumor], cols, rows, fill=None, edgecolor="red") # tumor
)
plt.imshow(slope3, cmap="jet")
plt.colorbar()
# plt.show()
plt.gca().add_patch(
plt.Rectangle([colNormal, rowNormal - rowOffset], cols, rows, fill=None, edgecolor="black") # normal
)
plt.gca().add_patch(
plt.Rectangle([colTumor, rowTumor - rowOffset], cols, rows, fill=None, edgecolor="red") # tumor
)
plt.subplot(1, 3, 1)
plt.hist(slope3[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), bins=100, facecolor=(0, 0, 0, 0.5)) # normal
plt.hist(slope3[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), bins=100, facecolor=(1, 0, 0, 0.5)) # tumor
plt.title("slope")
plt.subplot(1, 3, 2)
plt.hist(intercept3[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), bins=100, facecolor=(0, 0, 0, 0.5)) # normal
plt.hist(intercept3[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), bins=100, facecolor=(1, 0, 0, 0.5)) # tumor
plt.title("intercept")
plt.subplot(1, 3, 3)
plt.hist((slope3 * 16 + intercept3)[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), bins=100, facecolor=(0, 0, 0, 0.5)) # normal
plt.hist((slope3 * 16 + intercept3)[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), bins=100, facecolor=(1, 0, 0, 0.5)) # tumor
plt.title("midband")
plt.subplot(1, 3, 1)
plt.boxplot([
slope3[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), # normal
slope3[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1) # tumor
], showfliers=False, labels=["normal", "tumor"])
plt.title("slope")
plt.subplot(1, 3, 2)
plt.boxplot([
intercept3[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), # normal
intercept3[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1) # tumor
], showfliers=False, labels=["normal", "tumor"])
plt.title("intercept")
plt.subplot(1, 3, 3)
plt.boxplot([
(slope3 * 16 + intercept3)[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), # normal
(slope3 * 16 + intercept3)[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1) # tumor
], showfliers=False, labels=["normal", "tumor"])
plt.title("midband")
# plt.xticks([1, 2, 3], ["slope", "intercept", "midband"])
plt.tight_layout()
import scipy.io
namespace = scipy.io.matlab.loadmat("mouse2.mat")
slope2 = namespace["slope"]
intercept2 = namespace["intercept"]
waveform0 = plt.imread("mouse2-processed.png")
waveform = np.dot(waveform0[..., :3], [0.299, 0.587, 0.114])[0: 644, 11: 1023]
waveform
plt.imshow(waveform, cmap="gray")
plt.colorbar()
# plt.show()
plt.imshow(waveform, cmap="jet")
plt.colorbar()
plt.xticks([])
plt.yticks([])
# plt.show()
%%time
def task(i):
i, waveform = i
import numpy as np
import scipy.fftpack
fs = 300e+6
c0 = 1500
x0 = c0 / fs
t0 = 1 / fs
n = 30
nn = 2 * n + 1
# f = list(range(5, 29))
f = list(range(0, n+1))
hammingWindow = np.hamming(nn)
slope3 = np.zeros((waveform.shape[0] - nn + 1, 1))
intercept3 = np.zeros((waveform.shape[0] - nn + 1, 1))
def slopeFilter(X):
fft_temp = np.abs(scipy.fftpack.fft(X * hammingWindow))
FT = np.nan_to_num(20 * np.log10(fft_temp / np.max(fft_temp)))
# return np.polyfit(f, FT[5: 29], 1)
return np.polyfit(f, FT[0: n+1], 1)
for j in range(0, waveform.shape[0] - nn):
slope3[j], intercept3[j] = slopeFilter(waveform[j: j + nn])
return slope3, intercept3
# results = []
# f = list(range(0, n+1))
# for i in range(0, waveform.shape[1]):
# results.append(lview.apply(task, i, waveform[:, i]))
# print(i)
results = lview.map(task, [(i, waveform[:, i]) for i in range(0, waveform.shape[1])])
results
slope2 = np.hstack([i[0] for i in results.result()])
slope2
intercept2 = np.hstack([i[1] for i in results.result()])
intercept2
plt.imshow(slope2, cmap="jet")
plt.colorbar()
plt.xticks([])
plt.yticks([])
# plt.show()
plt.imshow(intercept2, cmap="jet")
plt.colorbar()
plt.xticks([])
plt.yticks([])
# plt.show()
plt.imshow(slope2 * 16 + intercept2, cmap="jet")
plt.colorbar()
plt.xticks([])
plt.yticks([])
# plt.show()
scipy.io.matlab.savemat("mouse2.mat", {
"slope": slope2,
"intercept": intercept2,
"midband": slope2 * 16 + intercept2
})
colNormal = 600
rowNormal = 200
colTumor = 620
rowTumor = 300
rows = 80
cols = 100
rowOffset = 30
plt.imshow(waveform, cmap="jet")
plt.colorbar()
# plt.show()
plt.gca().add_patch(
plt.Rectangle([colNormal, rowNormal], cols, rows, fill=None, edgecolor="black") # normal
)
plt.gca().add_patch(
plt.Rectangle([colTumor, rowTumor], cols, rows, fill=None, edgecolor="red") # tumor
)
plt.xticks([])
plt.yticks([])
plt.imshow(slope2 * 16 + intercept2, cmap="jet")
plt.colorbar()
# plt.show()
plt.gca().add_patch(
plt.Rectangle([colNormal, rowNormal - rowOffset], cols, rows, fill=None, edgecolor="black") # normal
)
plt.gca().add_patch(
plt.Rectangle([colTumor, rowTumor - rowOffset], cols, rows, fill=None, edgecolor="red") # tumor
)
plt.xticks([])
plt.yticks([])
plt.scatter(slope2[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), intercept2[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), c=(0, 0, 0, 0.2), marker=".")
plt.scatter(slope2[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), intercept2[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), c=(1, 0, 0, 0.2), marker=".")
plt.subplot(1, 3, 1)
plt.hist(slope2[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), bins=100, facecolor=(0, 0, 0, 0.5)) # normal
plt.hist(slope2[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), bins=100, facecolor=(1, 0, 0, 0.5)) # tumor
plt.title("slope")
plt.subplot(1, 3, 2)
plt.hist(intercept2[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), bins=100, facecolor=(0, 0, 0, 0.5)) # normal
plt.hist(intercept2[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), bins=100, facecolor=(1, 0, 0, 0.5)) # tumor
plt.title("intercept")
plt.subplot(1, 3, 3)
plt.hist((slope2 * 16 + intercept2)[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), bins=100, facecolor=(0, 0, 0, 0.5)) # normal
plt.hist((slope2 * 16 + intercept2)[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), bins=100, facecolor=(1, 0, 0, 0.5)) # tumor
plt.title("midband")
plt.subplot(1, 3, 1)
plt.boxplot([
slope2[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), # normal
slope2[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1) # tumor
], showfliers=False, labels=["normal", "tumor"])
plt.title("slope")
plt.subplot(1, 3, 2)
plt.boxplot([
intercept2[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), # normal
intercept2[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1) # tumor
], showfliers=False, labels=["normal", "tumor"])
plt.title("intercept")
plt.subplot(1, 3, 3)
plt.boxplot([
(slope2 * 16 + intercept2)[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), # normal
(slope2 * 16 + intercept2)[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1) # tumor
], showfliers=False, labels=["normal", "tumor"])
plt.title("midband")
# plt.xticks([1, 2, 3], ["slope", "intercept", "midband"])
plt.tight_layout()
plt.imshow((slope2 < -0.8) & (intercept2 < -11.25) & ((slope2 * 16 + intercept2) < -25), cmap="gray")
plt.colorbar()
plt.xticks([])
plt.yticks([])
import scipy.io
namespace = scipy.io.matlab.loadmat("Ivus_data.mat")
namespace
# waveform0 = namespace["Rabbit_aorta_80M_Sep8_pa_RF_01"]
waveform0 = namespace["Ivus_data"]
# waveform = np.dot(waveform0[..., :3], [0.299, 0.587, 0.114])[0: 644, 11: 1023]
# waveform = np.nan_to_num(np.log(np.abs(waveform0)))
waveform = waveform0
waveform
waveform.max(), waveform.min()
waveform.shape
plt.imshow(waveform, cmap="jet")
plt.colorbar()
# plt.show()
plt.imshow(waveform, cmap="gray")
plt.colorbar()
# plt.show()
%%time
def task(i):
i, waveform = i
import numpy as np
import scipy.fftpack
fs = 300e+6
c0 = 1500
x0 = c0 / fs
t0 = 1 / fs
n = 30
nn = 2 * n + 1
f = list(range(5, 29))
hammingWindow = np.hamming(nn)
slope3 = np.zeros((waveform.shape[0] - nn + 1, 1))
intercept3 = np.zeros((waveform.shape[0] - nn + 1, 1))
def slopeFilter(X):
fft_temp = np.abs(scipy.fftpack.fft(X * hammingWindow))
FT = np.nan_to_num(20 * np.log10(fft_temp / np.max(fft_temp)))
return np.polyfit(f, FT[5: 29], 1)
for j in range(0, waveform.shape[0] - nn):
slope3[j], intercept3[j] = slopeFilter(waveform[j: j + nn])
return slope3, intercept3
# results = []
# f = list(range(0, n+1))
# for i in range(0, waveform.shape[1]):
# results.append(lview.apply(task, i, waveform[:, i]))
# print(i)
results = lview.map(task, [(i, waveform[:, i]) for i in range(0, waveform.shape[1])])
results
slope4 = np.hstack([i[0] for i in results.result()])
slope4
intercept4 = np.hstack([i[1] for i in results.result()])
intercept4
plt.imshow(slope4, cmap="jet")
plt.colorbar()
# plt.show()
plt.imshow(intercept4, cmap="jet")
plt.colorbar()
# plt.show()
plt.imshow(slope4 * 16 + intercept4, cmap="jet")
plt.colorbar()
# plt.show()
scipy.io.matlab.savemat("rabbit.mat", {
"slope": slope4,
"intercept": intercept4,
"midband": slope4 * 16 + intercept4
})
import scipy.io
namespace = scipy.io.matlab.loadmat("rabbit_aorta_80MHz.mat")
namespace
# waveform0 = namespace["Rabbit_aorta_80M_Sep8_pa_RF_01"]
waveform0 = namespace["image_log"]
# waveform = np.dot(waveform0[..., :3], [0.299, 0.587, 0.114])[0: 644, 11: 1023]
# waveform = np.nan_to_num(np.log(np.abs(waveform0)))
waveform = waveform0
waveform
waveform.max(), waveform.min()
plt.imshow(waveform, cmap="gray")
plt.colorbar()
# plt.show()
plt.imshow(waveform, cmap="jet")
plt.colorbar()
# plt.show()
slope5 = np.fromfile("rabbit-slope.numpy").reshape(-1, 1000)
intercept5 = np.fromfile("rabbit-intercept.numpy").reshape(-1, 1000)
midband5 = slope5 * 16 + intercept5
slope5 = np.hstack([i[0] for i in results.result()])
slope5
slope5.tofile("rabbit-slope.numpy")
intercept5 = np.hstack([i[1] for i in results.result()])
intercept5
intercept5.tofile("rabbit-intercept.numpy")
plt.imshow(slope5, cmap="jet")
plt.colorbar()
# plt.xticks([])
# plt.yticks([])
# plt.show()
plt.imshow(intercept5, cmap="jet")
plt.colorbar()
# plt.xticks([])
# plt.yticks([])
# plt.show()
plt.imshow(slope5 * 16 + intercept5, cmap="jet")
plt.colorbar()
# plt.xticks([])
# plt.yticks([])
# plt.show()
colNormal = 100
rowNormal = 500
colTumor = 100
rowTumor = 200
rows = 80
cols = 100
rowOffset = 30
plt.imshow(waveform, cmap="jet")
plt.colorbar()
# plt.show()
plt.gca().add_patch(
plt.Rectangle([colNormal, rowNormal], cols, rows, fill=None, edgecolor="black") # normal
)
plt.gca().add_patch(
plt.Rectangle([colTumor, rowTumor], cols, rows, fill=None, edgecolor="red") # tumor
)
# plt.xticks([])
# plt.yticks([])
plt.imshow(slope5, cmap="jet")
plt.colorbar()
# plt.show()
plt.gca().add_patch(
plt.Rectangle([colNormal, rowNormal - rowOffset], cols, rows, fill=None, edgecolor="black") # normal
)
plt.gca().add_patch(
plt.Rectangle([colTumor, rowTumor - rowOffset], cols, rows, fill=None, edgecolor="red") # tumor
)
# plt.xticks([])
# plt.yticks([])
plt.imshow(intercept5, cmap="jet")
plt.colorbar()
# plt.show()
plt.gca().add_patch(
plt.Rectangle([colNormal, rowNormal - rowOffset], cols, rows, fill=None, edgecolor="black") # normal
)
plt.gca().add_patch(
plt.Rectangle([colTumor, rowTumor - rowOffset], cols, rows, fill=None, edgecolor="red") # tumor
)
# plt.xticks([])
# plt.yticks([])
plt.imshow(slope5 * 16 + intercept5, cmap="jet")
plt.colorbar()
# plt.show()
plt.gca().add_patch(
plt.Rectangle([colNormal, rowNormal - rowOffset], cols, rows, fill=None, edgecolor="black") # normal
)
plt.gca().add_patch(
plt.Rectangle([colTumor, rowTumor - rowOffset], cols, rows, fill=None, edgecolor="red") # tumor
)
# plt.xticks([])
# plt.yticks([])
plt.scatter(slope5[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), intercept5[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), c=(0, 0, 0, 0.2), marker=".")
plt.scatter(slope5[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), intercept5[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), c=(1, 0, 0, 0.2), marker=".")
plt.xlabel("slope")
plt.ylabel("intercept")
plt.subplot(1, 3, 1)
plt.hist(slope5[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), bins=100, facecolor=(0, 0, 0, 0.5)) # normal
plt.hist(slope5[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), bins=100, facecolor=(1, 0, 0, 0.5)) # tumor
plt.title("slope")
plt.subplot(1, 3, 2)
plt.hist(intercept5[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), bins=100, facecolor=(0, 0, 0, 0.5)) # normal
plt.hist(intercept5[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), bins=100, facecolor=(1, 0, 0, 0.5)) # tumor
plt.title("intercept")
plt.subplot(1, 3, 3)
plt.hist((slope5 * 16 + intercept5)[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), bins=100, facecolor=(0, 0, 0, 0.5)) # normal
plt.hist((slope5 * 16 + intercept5)[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1), bins=100, facecolor=(1, 0, 0, 0.5)) # tumor
plt.title("midband")
plt.subplot(1, 3, 1)
plt.boxplot([
slope5[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), # normal
slope5[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1) # tumor
], showfliers=False, labels=["normal", "tumor"])
plt.title("slope")
plt.subplot(1, 3, 2)
plt.boxplot([
intercept5[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), # normal
intercept5[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1) # tumor
], showfliers=False, labels=["normal", "tumor"])
plt.title("intercept")
plt.subplot(1, 3, 3)
plt.boxplot([
(slope5 * 16 + intercept5)[rowNormal - rowOffset: rowNormal + rows - rowOffset, colNormal: colNormal + cols].reshape(-1), # normal
(slope5 * 16 + intercept5)[rowTumor - rowOffset: rowTumor + rows - rowOffset, colTumor: colTumor + cols].reshape(-1) # tumor
], showfliers=False, labels=["normal", "tumor"])
plt.title("midband")
# plt.xticks([1, 2, 3], ["slope", "intercept", "midband"])
plt.tight_layout()